      SUBROUTINE CONVECT(OM,YM,PLOGADJ,TETA,UPI,VPI,PL,TECON,PCON,QPI)

!     Legacy Mars GCM v24
!     Mars Climate Modeling Center
!     NASA Ames Research Center

!
!  PURPOSE
!      CONVECT CHECKS FOR (SUPERADIABATIC) POTENTIAL TEMPERATURE
!      INSTABILITIES BETWEEN ADJACENT LAYERS FOR A GIVEN 'PI' GRID
!      POINT. ANY INSTABILITIES ARE RESOLVED BY THE PROCESS OF
!      ATMOSPHERIC CONVECTION. THE EXTENT AND THE POTENTIAL TEMPERATURE
!      OF THE CONVECTIVE ZONE ARE ALSO DETERMINED. (THE CONVECTIVE
!      ZONE STARTS AT THE SURFACE AND INCLUDES ALL LAYERS THAT EXCHANGE
!      HEAT WITH THE SURFACE THROUGH CONVECTION.)
!
!  AUTHOR
!      STEVE POHORSKY    INFORMATICS     TASK 605    OCT 82
!  FOR
!      JIM POLLACK (SEE NOTE OF 7-7-82. NOTE WAS WRITTEN BY BOB HABERLE)
!  ENVIRONMENT
!      Cray-2            UNICOS 3.0      FORTRAN
!  REVISION HISTORY
!      JR SCHAEFFER     TASK 904           DEC 86
!      INCORPORATE MOMENTUM EXCHANGE   (SEE NOTE OF 9/15/86).
!      Re-written September 1993.
!      Removed references to all commons and created an argument list
!      to pass all variables.
!  INPUT PARAMETERS
!      PL(K) ARRAY          - PRESSURE AT EACH LEVEL FOR THE 'PI' POINT.
!      TETA(K) ARRAY        - POTENTIAL TEMPERATURE AT EACH LEVEL.
!      OM(K) ARRAY          - CONVERSION FACTOR FROM TETA TO TL AT EACH
!                             LEVEL FOR THE 'PI' POINT.
!      YM(K) ARRAY          - THE MASS PER UNIT AREA OF EACH LAYER.
!                              (K = 2 TO NLEVSM1 BY 2S.)
!      PLOGADJ(K) ARRAY     - DIFFERENCE IN LOG OF PRESSURE BETWEEN
!                             ADJACENT LEVELS. PLOGADJ(K) = LN(PL(K+1))-
!                             LN(PL(K)). (K = 3 TO NLEVELS - 1)
!
!  OUTPUT PARAMETERS
!      TETA(K) ARRAY     - AS IN INPUT, BUT ADJUSTED FOR CONVECTION.
!      PCON              - PRESSURE AT THE TOP OF THE CONVECTIVE ZONE.
!      TECON             - POTENTIAL TEMPERATURE OF THE CONVECTIVE ZONE.
!
!  CALLED BY
!      COMP3
!
      use grid_h
      use defines_h
      implicit none

!######################################################################

!     IN CONVECT, THE TEMPERATURE AT THE MIDPOINT OF A LAYER IS USED
!     AS AN APPROXIMATION FOR THE TEMPERATURE OF THE WHOLE LAYER.
!     FOR THE ATMOSPHERE TO BE STABLE AT A POINT, THE POTENTIAL
!     TEMPERATURE MUST BE NON-DECREASING AS A FUNCTION OF ALTITUDE
!     (NON-INCREASING AS A FUNCTION OF THE LEVEL INDEX K FOR EVEN K.)
!     IF THE ATMOSPHERE IS UNSTABLE, CONVECT ADJUSTS THE POTENTIAL
!     TEMPERATURE VALUES UNTIL THE ATMOSPHERE IS STABLE. WHEN CONVECT
!     FINDS A SET OF ADJACENT LAYERS THAT LACKS STABILITY, IT CHANGES
!     THE POTENTIAL TEMPERATURE FOR ALL LAYERS IN THE SET TO AN
!     ENERGY-WEIGHTED AVERAGE VALUE IN A WAY THAT CONSERVES THE TOTAL
!     HEAT ENERGY OF THE LAYERS IN THE SET. ANY ADJUSTMENTS MADE THAT
!     INVOLVE LEVELS 4 OR 6 (THE TOP TWO LAYERS OF THE TROPOSPHERE)
!     INCLUDE AN ADJUSTMENT TO THE POTENTIAL TEMPERATURE OF THE
!     STRATOSPHERE ( TETA(3) ) TO KEEP TETA(3) A LINEAR EXTRAPOLATION
!     IN LOG OF PRESSURE FROM TETA(4) AND TETA(6). AS SPECIFIED IN THE
!     NOTE OF 7-7-82, THERE ARE THREE TYPES OF ADJUSTMENTS (TYPES A, B,
!     AND C) WHICH ARE USED DEPENDING ON THE INVOLVEMENT OF THE
!     STRATOSPHERE.  THE ALGORITHM USED IN CONVECT RESEMBLES
!     MATHEMATICAL INDUCTION.  START WITH SET OF THE TOP TWO LAYERS
!     OF THE TROPOSPHERE.  MAKE ANY ADJUSTMENTS REQUIRED TO PRODUCE
!     STABLILITY FOR THE SET.  ADD NEXT LAYER OF THE TROPOSPHERE TO
!     THE SET.  MAKE ANY ADJUSTMENTS REQUIRED TO PRODUCE STABLILITY
!     FOR THE SET.  CONTINUE ADDING LAYERS AND STABLIZING RESULTANT
!     SET UNTIL THE SET  CONTAINS ALL OF THE TROPOSPHERE.
!
!     IF THE BOTTOM LAYER GETS ADJUSTED IN THIS PROCESS, THEN THERE
!     WILL BE A (VERTICAL) INTERVAL FROM THE SURFACE OF THE PLANET TO
!     WHAT IS CALLED 'THE TOP OF THE CONVECTIVE ZONE' OVER WHICH THE
!     POTENTIAL TEMPERATURES ARE CONSTANT. IN THIS CASE, THIS INTERVAL
!     IS THE CONVECTIVE ZONE. (SEE BELOW IN CODE FOR THE OTHER CASE.)
!     (IF TWO ADJACENT LAYERS HAVE THE SAME POTENTIAL TEMPERATURE, THE
!     ADJUSTMENT MADE IS ACADEMIC AND ONLY SERVES TO ADVANCE THE
!     ALGORITHM IN DETERMINING THE EXTENT OF THE CONVECTIVE ZONE. THUS
!     THE BORDERLINE STABLE CASE IS HANDLED AS UNSTABLE.)
!
!#########

      real*8  :: WTC(L_LEVELS)

      REAL*8  :: OM(L_LEVELS), YM(L_LEVELS), TETA(L_LEVELS)
      real*8  :: PL(L_LEVELS)
      REAL*8  :: UPI(L_LEVELS), VPI(L_LEVELS), PLOGADJ(L_LEVELS)
! Include instantaneous tracer mixing by convection
      REAL*8  :: QPI(L_LEVELS,NTRACE), QCON(NTRACE)

      integer :: k, m, kk, kj, mkj, km2, kp3
      real*8  :: betacon, sum, sumuv, tecon, pcexp, pcon

!#=====================================================================

!     CALCULATE SOME FREQUENTLY USED EXPRESSIONS.

!     WTC is a weighting factor for convection.  WTC * TETA is
!     proportional to the heat energy of the layer.

      WTC(3) = OM(3) * YM(2)

      DO 100  K = 4, L_LEVELM1, 2
        WTC(K) = OM(K)*YM(K)
  100 CONTINUE

      BETACON = - PLOGADJ(3)/(PLOGADJ(4)+PLOGADJ(5))

!     EACH TIME THROUGH THE FOLLOWING LOOP ADDS ANOTHER LAYER TO THE
!     SET AND THEN MAKES THE ATMOSPHERE STABLE FOR THAT SET.

!     IF LEVELS K AND K+2 ARE STABLE, THE WHOLE SET OF LAYERS IS
!     STABLE. THIS IS BECAUSE IF K = 4, THERE ARE NO OTHER LAYERS, AND
!     IF K > 4, THE PREVIOUS TIME THROUGH THIS LOOP LEFT TETA(K)
!     NON-DECREASING IN ALTITUDE FOR LEVELS 4 TO K.

      DO 1000 K=4,L_LEVELM3,2

!       If stable, no adjustments are needed so by-pass the next
!       section

        IF (TETA(K).GT.TETA(K+2))   GOTO 900

!       INSTABILITY FOUND.  MAKE ADJUSTMENTS REQUIRED TO MAKE
!       TETA(K) NON-DECREASING IN ALTITUDE FROM LEVEL K+2 TO THE
!       TROPOPAUSE.

!       Set up for the next loop.

        SUM   = WTC(K+2)
        SUMUV = YM(K+2)
! INSERT TRACER MIXING IN VERTICAL (instantaneous)
        DO M = 1,NTRACE
          QCON(M) = QPI(K+2,M)
        END DO

        TECON = TETA(K+2)
        KP3   = K+3

!       WORK UPWARDS FROM LEVEL K+2 TO THE TROPOPAUSE INCLUDING ONE MORE
!       LAYER EACH TIME THROUGH THIS LOOP. THE FIRST TIME THROUGH THIS
!       LOOP, KJ = K AND TECON = TETA(K+2); SO WE ARE ADJUSTING THE
!       INSTABILITY WE FOUND ABOVE. THUS THE FIRST TIME THROUGH THIS
!       LOOP WE  GOTO 500 . AFTER THIS TETA(K) AND TETA(K-2) MAY NO
!       LONGER BE STABLE, SO WE KEEP WORKING UNTIL WE REACH THE
!       TROPOPAUSE OR UNTIL WE PERFORM AN ADJUSTMENT THAT DOES NOT CAUSE
!       ANY NEW INSTABILITY.
!
!       SINCE EACH TIME THROUGH THIS LOOP WE INCLUDE ONE MORE LAYER, THE
!       LOOP MUST HANDLE TWO CASES, THE CASE OF NO MORE INSTABILITIES
!       FOUND, AND THE CASE OF NEW INSTABILITIES FOUND.

        KM2   = K-2
        DO 700  MKJ = 2, KM2, 2
          KJ = K-(MKJ-2)

          IF (TETA(KJ).LE.TECON)   GOTO 500

!         CASE OF NO MORE INSTABILITIES FOUND. ONLY NEED TO CHECK FOR
!         PCON ADJUSTMENT. THEN EXIT THE MKJ LOOP.

          IF (TETA(KJ+1).GE.TECON)   GOTO 1000
          IF (K+2.NE.L_LEVELM1)      GOTO 1000

!         CHECK TETA AT LAYER BOUNDARY (INSTEAD OF LAYER MIDPOINT) AT
!         TOP OF CONVECTIVE ZONE. ADJUST PCON IF TETA NOT
!         NONDECREASING HERE.

          PCEXP = 2.0*(TECON-TETA(KJ+1))/(TETA(KJ)-TETA(KJ+1))
          PCON  = PL(KJ+1)*(PL(KJ)/PL(KJ+1))**PCEXP

          IF (PCON.LT.PL(KJ-1)) THEN
            PCON = PL(KJ-1)
          ENDIF

          GOTO 1000

500       CONTINUE

!         CASE OF NEW INSTABILITIES FOUND. THREE SUBCASES
!         ACCORDING TO KJ.

!         This 'if' statement prevents roundoff

          IF (TETA(KJ).NE.TECON)  THEN
            TECON = (TETA(KJ)*WTC(KJ)+TECON*SUM)/(WTC(KJ)+SUM)
          ENDIF

! TRACER MIXING
          DO M = 1, NTRACE
            IF(QPI(KJ,M).NE.QCON(M)) THEN
              QCON(M) = (QPI(KJ,M)*YM(KJ)+QCON(M)*SUMUV)/(YM(KJ)+SUMUV)
            ENDIF
          END DO

!         Adjust the layers from level KJ to level K.

          DO 600 KK = KJ, KP3
            TETA(KK) = TECON
! TRACER MIXING
            DO M = 1,NTRACE
              QPI(KK,M)  = QCON(M)
            END DO

  600     CONTINUE

          PCON  = PL(KJ-1)
          SUM   = SUM+WTC(KJ)
          SUMUV = SUMUV+YM(KJ)

  700   CONTINUE

!       Stop working if the tropopause is reached.

        GOTO 1000

!       CASE OF NO ADJUSTMENT REQUIRED WHEN LEVEL K+2 ADDED TO
!       THE SET. IF LEVEL K+2 IS THE MIDPOINT OF THE BOTTOM LAYER,
!       CHECK FOR A 'SHALLOW' CONVECTIVE LAYER. OTHERWISE REPEAT
!       THE K LOOP FOR THE NEXT LAYERR.

900     CONTINUE

        IF (K+2.NE.L_LEVELM1)   GOTO 1000

!       CHECK TETA AT LAYER BOUNDARY (INSTEAD OF LAYER MIDPOINT)
!       AT TOP OF CONVECTIVE ZONE. ADJUST PCON IF TETA NOT
!       NONDECREASING HERE.

        IF (TETA(L_LEVELM2).GT.TETA(L_LEVELM1))   GOTO 950

!       Shallow convective layer present.

        TECON = TETA(L_LEVELM1)
        PCEXP = 2.0*(TECON-TETA(L_LEVELM2))/
     *              (TETA(L_LEVELM3)-TETA(L_LEVELM2))
        PCON  = PL(L_LEVELM2)*(PL(L_LEVELM3)/PL(L_LEVELM2))**PCEXP

        IF (PCON.LT.PL(L_LEVELM4)) THEN
          PCON = PL(L_LEVELM4)
        ENDIF

        GOTO 1000

!       SHALLOW CONVECTIVE LAYER NOT PRESENT. WHOLE ATMOSPHERE
!       STABLE. NO CONVECTION FOUND.

!       Assume a little convection near the ground.

  950   CONTINUE

        TECON = 0.5*(TETA(L_LEVELS)+TETA(L_LEVELM1))
        PCON  = PL(L_LEVELM1)

 1000 CONTINUE

      RETURN
      END
